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Abstract 

Adjoint sensitivity calculation of stress, buckling and 
displacement constraints may be much less expensive 
than direct sensitivity calculation when the number of 
load cases is large. Adjoint stress and displacement 
sensitivities are available in the literature. Expressions 
for local buckling sensitivity of isotropic plate elements 
are derived in this study. Computational efficiency of 
the adjoint method is sensitive to the number of 
constraints and, therefore, the method benefits from 
constraint lumping. A continuum version of the 
Kreisselmeier-Steinhauser (KS) function is chosen to 
lump constraints. The adjoint and direct methods are 
compared for three examples: a truss structure, a simple 
HSCT wing model, and a large HSCT model. These 
sensitivity derivatives are then used in optimization. 

Introduction 

Sensitivity of structural response quantities, such as 
functions of stress or displacement components, to 
design parameters is useful in structural applications of 
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optimization, system identification, and probabilistic 
analysis. Finite difference derivatives are commonly 
used, but they are expensive computationally and prone 
to errors. Consequently, there is wide interest in 
analytical sensitivity methods, commonly classified as 
either direct or adjoint methods (e.g., Adelman and 
Haftka, 1989). 

The widely used direct method is obtained by direct 
differentiation of the equations of equilibrium, which 
are used to compute the structural responses or static 
solutions. The number of static solutions required in 
obtaining the derivatives is equal to the product of the 
number of design variables and the number of load 
cases. The adjoint method is obtained by differentiation 
of constraint functionals. The number of static solutions 
required to obtain the derivatives is equal to the number 
of constraints of interest. For problems with multiple 
load cases or multiple structural configurations, such as 
damaged versions of a structure, the adjoint method can 
be more efficient than the direct method. This 
computational advantage is enhanced for derivatives of 
stress functions because these derivatives can be 
calculated directly rather than obtained from derivatives 
of the displacements. Reduction of the number of 
constraints through constraint deletion and/or lumping 
enhances the efficiency of the adjoint method. 

The computational efficiency of the adjoint method is 
offset by greater implementation effort compared to the 
direct method, especially when based on differentiation 
of discretized structural models. Akgiin et al. (1998) 
showed that this implementation penalty is alleviated 
for stress constraints when the adjoint method is 
implemented based on the continuum equations. In this 
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case, the adjoint load is implemented as initial strains 
imposed in the elements of interest. 

Initial strains for the sensitivity of von Mises 
constraints with and without constraint lumping were 
derived in a previous work (Akgiin et ah, 1998). 
Implementation of the adjoint method was 
demonstrated for truss and plane-stress elements. A 
continuum version of the Kreisselmeier-Steinhauser 
(KS) functional that allowed lumping of the constraints 
with the adjoint method was investigated via a twenty- 
five-bar truss problem. The present study is a 
continuation of that work. 

The objectives of this study are to derive the adjoint 
initial strains for local buckling sensitivity computation 
and to demonstrate the efficiency of the adjoint method 
under multiple load cases with constraint lumping. It is 
shown that buckling constraints satisfy a homogeneity 
property that simplifies buckling sensitivity 
computation by obviating the need to compute partial 
derivatives of the buckling constraint. The example 
problems used are the structural optimization of a truss 
structure, a high-speed civil transport (HSCT) wing, 
and a much larger HSCT model. 


load cases (i.e., the number of vectors f) is n f , Eq. (2) 
has to be solved n p n f times. 

The adjoint method used here is based on a continuum 
equation and requires the imposition of an adjoint load 
in the form of an initial strain distribution. When the 
method is applied to a discrete model, an equation of 
the form 

T Cl n Cl 

Ku =f (3) 

has to be solved n c times where n c is the number of 
constraints of interest and, therefore, the number of 

adjoint load vectors f , u is the adjoint displacement, 
and f “ is the force equivalent to the adjoint load. For 
stress constraints, f is derived from an applied initial 
strain. A finite element code that has initial strain 
capability can easily generate the adjoint force, f" . The 
use of the adjoint displacement u” to calculate 
sensitivities is described subsequently 

Buckling sensitivity formulation will be carried out 
next. A perfect isotropic rectangular flat plate of 
thickness t and sides a and b may buckle under biaxial 
loading if the maximum value of 


Formulation of Buckling Sensitivity 

The discretized equations of equilibrium may be written 
in terms of the stiffness matrix K, the force vector f, 
and the displacement vector u as 



(4) 


Ku = f (1) 

Differentiation of this equation with respect to a design 
variable p gives 

Ku, = f' = -Ku ( 2 ) 


exceeds unity, where an applied stress resultant N t is 
positive when tensile, and 


A = 7T 2 t 3 H\ 


m 


2 A 


n 

2 TT 

a b 


with 


(5) 


where subscript p denotes derivative with respect to the 

p 

design variable, and f is called the pseudo load. The 
direct method solves for the displacement derivative 
from Eq. (2) and, if necessary, computes stress 
constraints from the displacement and its derivative. If 
the number of design variables is n p and the number of 


H = El 12(1 -V 2 ) (6) 

An isotropic infinite strip of width b having two simply 
supported sides, on the other hand, may buckle under 
pure in-plane shear if 
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i _ K > 2 

A, 52.7 t 3 * H 


(7) 


G - = iJk + *.«- (N .- +N .. £ - ) ]‘ M 


(13) 


exceeds unity (Whitney, 1987, p. 121). When Eq. (7) is 
applied to a rectangular flat plate, b is taken to be the 
shorter edge. Buckling under combined loading may 
occur if the value of the interaction equation 


1 _ 1 1 
A ~ + 

comb n s 


( 8 ) 


where a dot indicates inner product, a subscript 
following a comma shows a partial derivative, and a 
subscript without a comma denotes a total derivative. 
If p is the property of another structural member, the 
first two terms inside the square brackets in Eq. (13) are 
zero. Taking p=t , 


N, =De ; N e =rD. (14) 


reaches unity (Lekhnitskii, 1968). A buckling 
constraint for an isotropic plate under combined loading 
may be therefore written as 

g(N,t) = gl (N,t) + g 2 (N,t)-l<0 (9) 

where N is the vector of stress resultants and 

, 1 

= (1 °) 

K 

If A„ is negative, g, is taken to be zero. 


= g 2 ( N 


In this formulation, it can be shown (Haftka and 

Giirdal, 1992, pp. 312-316) that the adjoint load f 
needed to compute sensitivities is obtained from an 
initial strain given by g N /A, regardless of what p is, so 
that the loading depends on the buckling equation but 
not on the details of the finite element. The initial strain 
for the buckling constraint is thus given by (assuming 

K > o) 




(52.7 H) 2 t 6 


(15) 


The sensitivity of a buckling constraint is derived more 
conveniently if the constraint is expressed as 

G = ^jg(N,0<M (11) 

where A is the plate surface area. That is, G is the value 
of the buckling constraint averaged over an area A. . 
The constitutive equation for the plate is 



X" 


£ x 

N = 

N y 

= tD 

£ y 




Jxy_ 


where D is the constitutive matrix. The sensitivity of G 
to the design variable p , where p may, for example, be 
the plate's own thickness or the property of another 
structural element, is obtained by differentiating Eq. 
(11) with respect to p to yield 


This loading can be easily applied, then, in any program 
that has initial strain loading capabilities. Under the 
application of the adjoint load, the resulting 
displacement field is called the adjoint displacement u“. 
If we assume that g is constant over A it can be shown 
that 

G t =g t = 8j +g' K - N,+f'-u fl (16) 

where f ' is the pseudo load which is also used in the 
direct method (Eq. 2). 

To simplify Eq. (16), a property of homogeneous 
functions is introduced next with a brief digression. If h 
is a homogeneous function of degree n in the arguments 
N and t. that is, if 

h(cN,ct) = c"h(N,t) (17) 

and since 
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( 23 ) 


N 

t 



(18) G{g(N,t))= \-e pg{m) dA 

J A 


it can be shown that 


h t +h N -N, = y/t (19) 

(Akgiin et al., 1998). In the present case, the constraint 
g in Eq. (9) is not homogeneous as a whole but its 
components are. For g p n =~ 2 and for g y n=~ 4. 
Equation (16) hence becomes 

G,=g,=-j(2g 1 +4g 2 ) + f p -u a (20) 

The homogeneity property (Eq. 19) hence does away 
with the need to implement and compute partial 
derivatives of the buckling constraint, and allows us to 
use instead the constraint value, which is already 
available. When p is a property of a structural element 
other than the plate whose sensitivity is being 
computed, the buckling sensitivity is simply given by 

G p =g p =f p -u a ( 21 ) 

It should be noted that the homogeneity is not limited to 
the constraints for isotropic plates; it is valid for 
composite laminates as well. 


Here the integral is taken over all the plate elements 
being lumped together, and A is the surface area of the 
corresponding plate. In the implementation for this 
study the integrand in Eq. (23) is constant within each 
plate. The derivative of the KS function is given as 

KS p = G p /pG (24) 


Steps similar to those leading to Eqs. (16) and (20) 
from Eq. (10) yield 




P~ 


(25) 


where m is the number of buckling constraints lumped 
together and (5 tip is a Kronecker delta. (5 Pp is unity if 

design variable p is the thickness of the ith plate and 
zero if design variable p is not the thickness of the ith 
plate. The adjoint displacement u“ is now due to an 
initial strain state, imposed simultaneously in each of 
the plate elements being lumped. The strain state to be 
imposed in the ith plate is given by 


pe 


pg, 


8 , N 


(26) 


Constraint Lumping 

The adjoint method requires computation of the adjoint 
displacements for every potentially active constraint. 
Reducing the number of constraints by lumping them 
into groups can increase the savings brought about by 
the method. The continuum version of the KS (Akgiin 
et al., 1998) function is used here to lump a set of 
buckling constraints for a number of flat plates into a 
single constraint and is defined here as 

KS(g{N,t)) = UnG(g(Nj)) ( 22 ) 

with p being a positive parameter and where 


where g N /A is given by the right hand side of Eq. (15), 
which is the strain state imposed in a single element in 
the absence of lumping. 

When multiple load cases are present, one may want to 
lump all constraints for the group of m members under 
all the load cases together into the same KS function. 
This is the approach we used with the direct method in 
the comparison studies described below. For such an 
approach, however, we could not find an initial strain 
state for sensitivity calculation with the adjoint method. 
Instead of lumping, approach adopted here for the 
adjoint method is to select the most critical load case 
for buckling and von Mises constraints separately for 
each group of elements and to lump each type of 
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constraints under the respective load case together. The 
constraints under the other load cases are ignored. The 
initial strain state given by Eq. (26) is then computed at 
the most critical load case. 


Applications 

Three examples of varying complexity are used to 
compare the efficiency of the adjoint and direct 
methods. The three examples are a 108-bar truss, a 
simple HSCT wing model, and a large HSCT model. 
Sensitivity derivatives are obtained using both methods 
and then used in optimization. The 108-bar truss 
example is used to compare the effect of lumping on 
the two methods. The efficiency of both methods is 
compared by varying the number of load cases for the 
two HSCT examples. 

Both the adjoint and direct methods are implemented in 
a commercial finite element program, EAL 
(Engineering Analysis Language, EISI-EAL, 1983), for 
plane-stress elements with von Mises stress constraints 
and buckling constraints and for truss elements. The 
optimization code MINOS (Murtagh and Saunder, 
1993) has been incorporated into EAL. The linear 
programming option of MINOS (based on the simplex 
method) is used in a sequential linear programming 
(SLP) optimization approach for this study. 

Examples 1 and 2 were run on a Pentium 166 MHz PC 
with a single processor with a Windows NT operating 
system at the University of Florida. Example 3 was run 
on a dual processor Pentium 400 MHz PC with a 
Windows NT operating system at the NASA Langley 
Research Center. For all three examples, a finite 
difference approach was used to compute the derivative 
of the stiffness matrix, design variables were 
constrained by move limits which were initially 50 
percent of the starting design variable values, and when 
lumping was used the value of the KS function 
parameter p was 50. 

The lumping strategy used with the adjoint method was 
to lump each type of constraints (i.e., von Mises and 
buckling types where applicable) in a design region 
under the most critical load case into a single KS 


function whereas the strategy with the direct method 
was to lump each type of constraints under all the load 
cases into a single KS function. The two lumping 
strategies were thus different even though the resulting 
optima were very close. No constraint deletion strategy 
was used with either method with or without lumping. 

Example 1: 108-bar truss 

A 45-node, 108-bar, 80-dof plane truss, shown in Fig. 
1, was divided into 32 design regions, with bars in each 
region having a common cross-section, thus giving 32 
design variables. Optimization runs were made with 
seven load cases using the direct and the adjoint 
methods, with and without lumping. When lumping 
was used, bar stress constraints in each region were 
lumped into a continuum KS function. Hence there 
were 32 KS functions in either method. When no 
lumping was used, there were 756 constraints. The 
direct method with and without lumping required the 
solution of Eq. (2) 224 times, with the stiffness matrix 
already factored, whereas the adjoint method required 
the solution of Eq. (3), of the same size, 108 and 32 
times without and with lumping, respectively. The 
reason for 108 solutions instead of 756 is that the 
adjoint load vector f" for a given bar for different load 
cases differs only by a scalar factor. 

The optimum weight and the CPU time using the 
adjoint method and the direct method are given in Table 
1. The CPU times are for the entire optimization run 
(i.e., 8 SLP cycles). The optimum weights in the table 
are the weights at the end of eight optimization cycles 
representing convergence to 10 significant digits. 

As can be seen from Table 1, the final weights show 
some dependence on lumping and the choice of 
method. This is attributed to optimization history and to 
the effect of the KS function, which applies a slightly 
more conservative constraint. It can also be seen that 
the lumping was more beneficial to the direct method 
than the adjoint method, even though one would expect 
the opposite. The reason for this anomaly is that the 
savings associated with the adjoint method is in the 
solution of the equations of equilibrium with fewer 
right hand sides. This solution increases superlinearly 
with problem size, and is very small for small 


5 

American Institute of Aeronautics and Astronautics 



problems. The additional cost due to implementation 
discussed earlier, is proportional to the number of 
elements in the problem, and is therefore more 
significant than the solution time for small problems. 

Table 1. Optimum weights of the 108-bar truss and 
total (8 cycles) CPU times with the direct and 
adjoint methods with and without KS lumping 
under seven load cases. 



Weight, lb 
(Time, sec) 

Direct 

No 

Lump 

36.523367 

(312) 

Lump 

37.359306 

(169) 

Adjoint 

No 

Lump 

37.022500 

(292) 

Lump 

37.353403 

(256) 


Example 2: Simple HSCT wing 

A simple high-speed civil transport (HSCT) wing 
model developed by Balabanov et al. (1996), whose 
planform view is shown in Fig. 2, is investigated here. 
The 3-D model has 193 nodes, 449 2-node bar 
elements, 383 triangular isotropic membrane elements 
and 129 isotropic shear panel elements with a total of 
533 degrees of freedom. Figure 2 shows the design 
regions selected on the top skin of the HSCT model, the 
bottom skin being the mirror image of those. Skin 
regions 1, 2, 10, 11, 12, and 13, which cover the trailing 
edge and strake of the wing, are governed by the 
minimum gauge-thickness constraint set at 0.0035 ft. 
Skin thickness within each design region is uniform. 
The minimum gauge area used for the bar elements is 
0.0040 ft 2 . Not all of the finite elements are designed. 
There are 40 design variables consisting of 26 
membrane thicknesses and 14 bar cross-sectional areas. 

A series of optimization runs with different numbers of 
load cases was made. There were five basic load cases 


which were due to pull-up maneuvers, climb, cruise and 
taxiing. Additional fifteen load cases were derived as a 
combination of the basic five cases. Stress constraints 
were based on von Mises equivalent stresses in the 
membrane elements and axial stresses in the bars. 

The implementation of the adjoint method for the 
sensitivity of lumped buckling constraints was first 
verified. The model was then optimized under one, five 
and twenty load cases in different runs subject to 
buckling and von Mises stress constraints for 
membrane elements and yield limits for bar elements. 
There were 66 KS functions with both methods 
regardless of the number of load cases. Equation (2) 
had to be solved 40 n f times with the direct method 
where n f is the number of load cases and Eq. (3) 66 
times with the adjoint method regardless of the number 
of load cases. The optimum weights and the CPU times 
are given in Table 2. The CPU times are for the entire 
optimization run (i.e., 8 SLP cycles). The optimum 
weights in the table are the weights at the end of eight 
optimization cycles representing convergence to 13 
significant digits. 

Table 2. Optimum weights of the simple HSCT 
model and total (8 cycles) CPU times with the two 
methods using KS lumping. 


No of 

Direct 

Adjoint 

Load 

W, lb 

W, lb 

Cases 

(t, sec) 

(t, sec) 

1 

119746.87495532 

119746.87495534 


(750) 

(1410) 

5 

149327.13548552 

149327.13548557 


(1410) 

(1560) 

20 

149327.13548552 

149327.13548557 


(2962) 

(1627) 


As can be seen from Table 2, the final weights agree to 
13 digits accuracy. Also, for this larger problem, the 
improved efficiency with number of load cases of this 
implementation of the adjoint method compared to the 
direct method is apparent. For five load cases the two 
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methods are comparable, and for 20 load cases the 
adjoint method is about 45% cheaper. 

Example 3: Large HSCT 

The third example is a half-symmetric model of an 
entire high-speed transport aircraft 1 shown in Fig. 3. 
This model was originally presented by Scotti (1995). 
Results were performed at NASA Langley Research 
Center. The present EAL adjoint runstreams cannot 
handle four-noded membrane elements and, therefore, 
each four-noded element in the original model’s design 
regions was replaced with two three-noded elements, 
which increased the total number of elements in the 
model to 16434. In addition, all composite honeycomb 
plate sections were replaced with isotropic membrane 
sections, which reduced the number of design variables 
in the model to 61. Both stress and buckling constraints 
are computed for the triangular elements in the design 
regions. 

There are seven basic load cases, six arising from 
supersonic maneuvers at +2.5g and -l.Og and the 
seventh represents a taxi condition. Additional fourteen 
load cases were generated from various combinations 
of the basic six maneuver load cases. For timing 
purposes, the model was run for four optimization 
cycles using both the adjoint and direct method. For 
each method, analyses were performed with one, seven, 
14, and 21 load cases. The CPU time required to 
complete an optimization cycle is plotted against the 
number of load cases in Fig. 4. The lines shown in the 
figure are best-fit through the average of the run times 
for the four cycles. The actual run times for individual 
cycles are also indicated on the figure. 

For a single load case, the direct method is twice as fast 
as the adjoint method. However, as the number of load 
cases increases, the adjoint method becomes more 
efficient. Because this example is larger than the prior 
two examples, the efficiency of the adjoint method as 


1 The structural model for this example has been 
supplied by the Boeing Company and the results are 
presented without absolute scales in this paper under 
the conditions of a NASA Langley Property Loan 
Agreement, Loan Control Number 1922931. 


the number of load cases increases is greater, and the 
break-even point now appears to be around three or 
four load cases. For 21 load cases, the problem 
solution with the adjoint method requires only about 
one-third of the time required by the direct method. 

In addition, the optimization was run to convergence 
using both methods with one load case and seven load 
cases. Figure 5 shows the convergence history for the 
normalized objective function. It is seen that the direct 
and adjoint methods converge to comparable weights. 
The values of the design variables were also similar. 

Discussion and Conclusions 

Sensitivity of normal and shear buckling constraints in 
membranes has been formulated and implemented 
using the adjoint method. Sensitivity of lumped 
constraints has also been implemented. The present 
formulation utilizes a characteristic of homogenous 
functions to avoid calculation of partial derivatives of 
the constraints and the shear resultants. Test cases 
including a 108-bar truss and two high speed civil 
transport (HSCT) models of different complexities 
were used to compare the direct and the adjoint 
methods of sensitivity computation. The two methods 
have been compared in terms of efficiency for various 
numbers of load cases in optimization runs performed 
with the EAL finite element software. 

The direct method was found to be more efficient for a 
small number of load cases. The adjoint method 
outperforms the direct method in runs with a large 
number of load cases. The advantage of the adjoint 
method results from requiring static solutions equal in 
number to the number of active constraints, while the 
direct method requires a number equal to the number of 
design variables times the number of load cases. The 
advantage of the adjoint method increases with problem 
size because the solution cost per load case increases 
superlinearly with problem size. It must be noted, 
however, that the implement EAL (Engineering 
Analysis Language) optimization code was originally 
developed for the direct method and the adjoint method 
has been recently added to the code; there is room for 
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improving the efficiency of the adjoint method as 
currently implemented. 

The implementation of the direct method in EAL has 
been aided with specialized processors developed 
especially for this method. The implementation of the 
adjoint method used standard EAL runstream 
commands, which may accrue higher overhead during 
execution. Thus, for a single load case the adjoint 
method required substantially more time than the direct 
method. 

No constraint deletion strategy was used with either the 
direct or the adjoint method. It may be worthwhile to 
investigate the effect of constraint deletion in 
conjunction with lumping or without lumping. 
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Figure 2. Simple HSCT wing model and the design regions on the top skin. 
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CP seconds/optimization cycle 



Figure 3. Half- symmetric finite-element model of the large high-speed civil transport. 



Figure 4. CPU timing as a function of load cases for the large high-speed civil transport. 
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Normalized objective function 



Figure 5. Optimization history for the large high-speed civil transport. 
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